#! /usr/bin/env python
"""
editor for NM64 data

Copyright (C) 2010-2013 Christian T. Steigies <steigies@physik.uni-kiel.de>

This module offers editor functions for NM64 data.
The functions are converted from the IZMIRAN Pascal module U_Editor.pas

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see <http://www.gnu.org/licenses/>.
"""

__author__ = "Christian T. Steigies <steigies@physik.uni-kiel.de>"
__license__ = "GPL License"

## --------------------------------------------------------------------------
__all__ = ["editor_simply_summa", "editor_median_summa", "editor_median", "editor"]

## --------------------------------------------------------------------------
def perturbation(nt, efp, k):
    """pertubation"""
    nt[k] = (1-efp/100)*nt[k]

    return (nt)

## --------------------------------------------------------------------------
def editor_simply_summa(kmax, Nt):
    """simply summa"""
    m = 0; NS = 0
    for i in range(0, kmax):
        if Nt[i] != 0:
            NS += Nt[i]
            m += 1

    if m != 0:
        NS = kmax*NS/m
        
    return (NS)

## --------------------------------------------------------------------------
def editor_median_summa(k, n0, e0):
    """median summa"""
    ik = 0
    sn = 0.0
    r = d = et = nt = [0.]*k
    tet = [0.]*3
    
    for i in range(0, k):
        if nt[i] > 0.0:
            ik += 1
            r[i] = nt[i] / e0[i] / n0[i]
            d[i] = r[i]
            et[i] = r[i]
        if ik < 3:
            sn = 0.0
            for i in range(0, k):
                et[i] = 0.0

        jk = ik / 2 +1 # Note: Ik DIV 2 + 1
        for j in range(0, jk):
            c1 = 0.0
            for i in range(0, k):
                if r[i] >= c1:
                    c1 = r[i]
                    ii = i
            if j < jk:
                c2 = c1
                r[ii] = 0.0
        tet[1] = c1
        if jk == ik:
            tet[1] = (c1+c2)/2.0
            for i in range(0, k):
                r[i] = abs(d[i]-tet[1])
            for j in range(0, jk):
                c1 = 0.0
                for i in range(0, k):
                    if r[i] >= c1:
                        c1 = r[i]
                        ii = i
                if j < jk:
                    c2 = c1
                    r[ii] = 0.0
            tet[2] = c1
            if jk == ik:
                tet[2] = (c1+c2)/2.0
                for i in range(0, k):
                    r[i] = d[i]
                for i in range(0, k):
                    et[i] = et[i] / tet[1]
                    sig = tet[2]
                    if sig < 0.015:
                        sig = 0.015
                        sig = 4*sig
                        m = 0
                        sn = 0
                for i in range(0, k):
                    if (abs(r[i]-tet[1]) < sig):
                        sn += nt[i]
                        m += 1
                    if (m < k/4): # k div 4
                        sn = 0
                    else:
                        sn = k*sn/m
                        
        
    return(nt, et, sn)

## --------------------------------------------------------------------------
def editor_median(k, n0, e0, nt):
    """median"""
    i = j = ik = jk = ii = 0
    c1 = c2 = s0 = sn = tet = 0.0
    # define all lists seperately, otherwise the same list will have different names!
    ne = [0.]*k
    et = [0.]*k
    r = [0.]*k
    
    for i in range(0, k):
        if nt[i] > 0.0:
            ik += 1 # number of counters with countrate >0
        s0 += n0[i]/e0[i] # sum, corrected by efficiency
        r[i] = nt[i] / e0[i] / n0[i] # ratio countrate versus reference value, corr by eff
        et[i] = r[i] # save efficiency for each counter

    if ik < 3: # if less than three counters, return 0, EXIT
        sn = 0.0
        for i in range(0, k):
            et[i] = 0.0
            #exit ???
            
    jk = ik/2 + 1 # Note: Ik DIV 2 + 1
    ## find median of ratio
    for j in range(0, jk):
        c1 = 0.0
        for i in range(0, k): # find max ratio and index
            if r[i] >= c1:
                c1 = r[i] 
                ii = i
        if (j+1) < jk: # for last value (jk) set ratio of max to 0.0
            c2 = c1
            r[ii] = 0.0 

    tet = c1 # ratio of max
    jk = (ik / 2) *2    
    if jk == ik: # even number of counters with countrate > 0.0
        tet = (c1+c2)/2.0 # tet = averate of ratio of max and second last counter
    
    sn = 0
    for i in range(0, k):
        et[i] /= tet
        ne[i] = nt[i] / et[i]
        sn += ne[i] # or sn = s0*tet

    return(nt, ne, et, sn)

## --------------------------------------------------------------------------

## --------------------------------------------------------------------------
#def editor(k, sf, ck, q, n0, e0, nt, ne, et, su, sd, se, ns, sn, sc, es, pd, ind):
def editor(k, n0, e0, nt):
    """SuperEditor
    algorithm described in Belov at el
    k      - number of channels
    n0[k]  - initial count rates (baseline)
    e0[k]  - initial efficiencies
    nt[k]  - current countrates 
    ne[k]  - corrected countrates
    et[k]  - calculated efficiency
    su[k]  - confidence interval: upper boundary
    sd[k]  - confidence interval: lower boundary
    se[k]  - standard deviation of efficiencies in percent
    ind[k] - channel condition: 1 ok, 0 failed
    ns     - sum of corrected countrate
    sn     - standard deviation of sum of corrected countrate
    sc     - theoretical standard deviation of corrected countrate
    es     - detector efficiency
    pd     - detector efficieny of original data
    """
    
    q = 3.
    ck = 1. #???
    sf = 1.
    ps = 1.
    from math import exp, log, sqrt
    
    print "n0\t", n0
    print "e0\t", e0
    print "nt\t", nt
    print ""

    ne = [0.]*k; et = [0.]*k; su = [0.]*k; sd = [0.]*k; se = [0.]*k; ind = [0]*k
    ns = sn = sc = es = pd = 0.
    j = km  = 0
    t = t1 = t2 = s2 = see = d = r = 0.
    
    # ------------
    def NUL():
        """NUK nt, is this global?"""
        for i in range(0, k):
            nt[i] = 0.0
#            et[i] = se[i] = ne[i] = sd[i] = su[i] = 0.0
#            ind[i] = 0
        return()
    # ------------

    lm = -1
    for m in [1, 2]: # two runs: estimation and final
        t0 = 0
        for i in range(0, k):
            ind[i] = 1
            if nt[i] < 1:
                ind[i] = 0 # defective channels
                t0 += n0[i]/e0[i]

        ## main loop
        REPEAT = True
        while(REPEAT):
            dm = 0
            for i in range(0, k):
                e1 = e2 = s1 = l = 0
                if nt[i] >= 1:
                    for j in range(0, k):
                        if not((ind[j]==0) or (j==i)):
                            print nt[i], n0[i], nt[j], n0[j]
                            print (1.*nt[i]/n0[i]), log(1.*nt[j]/n0[j])
                            r  = log(1.*nt[i]/n0[i]) - log(1.*nt[j]/n0[j])
                            l += 1; e1 += r; e2 += r*r; s1 += 1/n0[j]
                            
                    if l <= 1:
                        NUL()
                        print "NUL, exit(1)"
                        exit(1) # not enough valid channels, exit
                        
                    et[i] = e1/l
                    see = 0
                    if l > 1: # always true?
                        see = (e2-e1*et[i])/l/(l-1)
                    et[i] = exp(et[i])
                    if m == 2:
                        se[i] = 100*et[i]*sqrt(t0/ns*ck/sf/n0[i] + see)
                    sd[i] = q*e0[i]*sqrt(ck/sf*(1/n0[i] + s1/l/l))
                    if m == 2:
                        sd[i] *= sqrt(t0/ns)
                        
                    if ind[i] > 0:
                        d = abs(log(et[i]) - log(e0[i]))
                        if d >= dm: # find extreme value of efficiency
                            lm = i
                            dm = d

## wiederherstellung der daten defekter kanaele
## zaehlratenkalkulation (NS) des monitors und desssen effizienz (es)
## die effizienz von defekten kanaelen erhaelt ein negatives vorzeichen

            print "n0\t", n0
            print "e0\t", e0
            print "nt\t", nt
            print "ind\t", ind
            print "ne\t", ne
            
            l = sum(ind)
            if lm < 0: # was lm = 0, 0 is a valid index, use -1?
                ns = sn = sc = es = pd = 0
                break #goto MetkaNul

            if (dm < sd[lm]) and (l>0):
                break # goto Metka
            else:
                ind[lm] = 0
                if (2*l-k-1) < 0:
                    NUL()
                    print "NUL(), EXIT"
                    
            REPEAT = ((2*l-k-1) < 0)
##Metka:
        if lm >= 0:
            if ind[lm] != 0: # == 1?
                l += 1
               
        ns = t = t2 = e1 = e2 = 0
        for i in range(0, k):
            su[i] = e0[i]+sd[i]
            sd[i] = e0[i]-sd[i]
            t += nt[i]
            if ind[i] != 0: # valid channel
                ne[i] = nt[i]/e0[i]
                r  = nt[i]/n0[i]
                e1 += r
                e2 += r*r
                t2 += nt[i]/e0[i]
            else: # invalid channel
                et[i] = -et[i] # mark with negative efficiency
                l = t1 = 0
                for j in range(0, k): # sum of all valid channels
                    if not((j==i) or (ind[j]==0)):
                        t1 += nt[j]/n0[j]/e0[j]
                        l += 1
                ne[i] = t1*n0[i]/l # replace invalid channel with average from valid
            ns += ne[i] # sum of all channels (including corrected)

# MetkaNul false
    if lm >= 0:
        es = t/ns
        sc = 100*ck/sqrt(sf*t2)
        sn = 0
        if l > 1:
            sn = 100*sqrt(t0/ns*abs(e2-e1*e1/l)/l/(l-1))
        pd = t2/ns

#MetkaNul:
    for i in range(0, k):
        if nt[i] == 0:
            et[i] = su[i] = sd[i] = se[i] = ne[i] = 0

    return(nt, ne, et, su, sd, se, ns, sn, sc, es, ps, ind)

## --------------------------------------------------------------------------
if __name__ == "__main__":
    print "This is a module to edit NM64 data."

    kmax = 6
    SF = 60. #???
    SF = 1.
    N0 = [597, 674, 575, 556, 592, 534] # baseline values
    Nok = [595, 677, 579, 546, 592, 525] # current values
    Nt = [595, 677, 579, 546, 533, 525] # current values, after 10% decrease in ch#5
    E0 = [1.0, 1.0, 1.0, 1.0, 1.0, 1.0] # efficiency
    pressure = 991.84
    multiplicity = 1.41
    nsigma = 3

    print "N0\t", N0, "\t", sum(N0)
    print "Nt ok\t", Nok, "\t", sum(Nok)
    print "Nt\t", Nt, "\t", sum(Nt)
    print

    print "Nin\t", Nt,
    NS = editor_simply_summa(kmax, Nt)
    print "\t", NS, "SumEditor"
    print

    nt = [0]*kmax
    for i in range(0, kmax):
        nt[i] = SF*Nt[i]

#    print "perturb", nt
#    nt = perturbation(nt, 10, 5) perturbation does not work!
#    print "perturb", nt
    (nt, ne, et, sn) = editor_median(kmax, N0, E0, nt)
    ns = int(sn / SF + 0.5)
    for i in range(0, kmax):
        nt[i] = int(nt[i]/SF + 0.5)
        ne[i] = int(ne[i]/SF + 0.5)
        # et[i] = (int(1000*et[i] + 0.5) / 1000.)

    print "Nin\t", nt
    print "Nout\t", ne,
    print "\t", ns, "medianEditor"
    print "Eff\t", et
    print

    (nt, ne, et, su, sd, se, ns, sn, sc, es, ps, ind) = editor(kmax, N0, E0, nt)
    print "Nin\t", nt
    print "Nout\t", ne,
    print "\t", ns, "SuperEditor"
    print "Eff\t", et
    print "+3s\t", su
    print "-3s\t", sd
    print


"""nmforwww.opt
D:\STATION\Kiel\PRIMARY
no D:\STATION\Kiel\DATA_NM
C:\INETPUB\FTPROOT\COSRAY!\FTP_NM
    19      0.00000000    0.10000000    0.00000000
    20    -35.44000000    0.11600000   -0.00004124
    21    -34.09000000    0.11658000   -0.00004388
    22      8.79000000   -0.00164000    0.00000046
    23      8.82000000   -0.00158000    0.00000041
    24      8.82000000   -0.00158000    0.00000041
     0      0.00000000    0.00000000    0.00000000
  30220  35780  36960  37362  37020  31690
  32550  36785  38000  37284  36826  32270
  32200  37003  37685  37910  37125  31030
 -30000  30000  30000  30000  30000  30000
  1.000  1.000  1.000  1.000  1.000  1.000
  1.000  1.000  1.000  1.000  1.000  1.000
  1.000  1.000  1.000  1.000  1.000  1.000
  1.000  1.000  1.000  1.000  1.000  1.000
1006.60   0.00   3.00
  -0.72   1.41   0.00
   0.00 100.00   1.00
   0.00   0.00   0.00
 1.0000
18.0000
 0.0000
 1.0508
 2.0000
 0.0000
 7.0000
 0.0000
 1.0000
 0.0000
 0.0000
24.0000

"""
